#include<iostream>
#include <xmmintrin.h> //SSE
#include <emmintrin.h> //SSE2
#include <pmmintrin.h> //SSE3
#include <tmmintrin.h> //SSSE3
#include <smmintrin.h> //SSE4.1
#include <nmmintrin.h> //SSSE4.2
#include <immintrin.h> //AVX、AVX2
#include<math.h>
#include<windows.h>
#include<omp.h>
using namespace std;
int N = 2000;
const int NUM_THREADS = 7;

// 生成不会出错的测试用例,事实是会出错
void m_reset(float** m, int n)
{
    for (int i = 0; i < n; i++)
    {
        // 下三角全部设置为0
        for (int j = 0; j < i; j++)
            m[i][j] = 0;
        //对角线设置为1.0
        m[i][i] = 1.0;
        //上三角设置为随机数
        for (int j = i + 1; j < n; j++)
            m[i][j] = (rand() % n);
    }

    //每一行加上比自己下标小的行
    for (int i = 0; i < n; i++)
    {
        int row1 = rand() % n;
        int row2 = rand() % n;
        int row3 = rand() % n;
        for (int j = 0; j < n; j++)
        {
            m[i][j] = m[i][j] + 2 * m[row1][j] + m[row2][j] + m[row3][j];
        }

    }

}

// 检查矩阵的正确性
void show(float** m, int n)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
            cout << m[i][j] << " ";
        cout << endl;
    }
}

// 平凡
void normal_solve(float** m, int n)
{
    #pragma omp parrallel if(parallel),num_threads(NUM_THREADS),private(i,j,k,tmp)
    for (int k = 0; k < n; k++)
    {
        //处理第k行
        //均除以对应行位于对角线上的元素
        #pragma omp single
        {
            __m256 vt = _mm256_set1_ps(m[k][k]);
            int j;
            //主线程做的事情，处理本行数据
            for (j = k + 1; j + 8 < n; j = j + 8)
            {
                //matrix[k][j] = matrix[k][j] / matrix[k][k];
                //非对齐模式
                __m256 va = _mm256_loadu_ps(&m[k][j]);
                va = _mm256_div_ps(va, vt);
                _mm256_storeu_ps(&m[k][j], va);
            }
            //处理剩余的
            if (j < n)
            {
                for (; j < n; j++)
                    m[k][j] = m[k][j] / m[k][k];
            }
            m[k][k] = 1.0;
        }

        #pragma omp for
        for (int i = k + 1; i < n; i++)
        {
            //串行原理的理解，m[k][k]变成了1，则第i行减去的次数为
            //m[i][k]/m[k][k]=m[i][k]
            int j;
            __m256 vaik = _mm256_set1_ps(m[i][k]);
            for (j = i + 1; j + 8 <n; j = j + 8)
            {
                __m256 vakj = _mm256_loadu_ps(&m[k][j]);
                __m256 vaij = _mm256_loadu_ps(&m[i][j]);
                __m256 vx = _mm256_mul_ps(vakj, vaik);
                vaij = _mm256_sub_ps(vaij, vx);
                _mm256_storeu_ps(&m[i][j], vaij);
            }
            //处理剩余的
            if (j < n)
            {
                for (; j < n; j++)
                    m[i][j] = m[i][j] - m[k][j] * m[i][k];
            }
            m[i][k] = 0.0;
        }
    }
    return;
}


int main()
{
    float** matrix_normal;
    matrix_normal = new float* [N];
    for (int i = 0; i < N; i++)
    {
        matrix_normal[i] = new float[N];
    }

    int size = 32;//数据集增长速度
    long long head1, tail1, freq1;//记录串行时间
    for (int n = size; n < N; n += size)
    {
        m_reset(matrix_normal, n);

        //串行
        QueryPerformanceFrequency((LARGE_INTEGER*)&freq1);
        QueryPerformanceCounter((LARGE_INTEGER*)&head1);
        QueryPerformanceCounter((LARGE_INTEGER*)&tail1);
        int count1 = 0;
        while ((tail1 - head1) < 0.1 * freq1) {
            normal_solve(matrix_normal, n);
            QueryPerformanceCounter((LARGE_INTEGER*)&tail1);
            count1++;
        }
        float time1 = (tail1 - head1) * 1000.0 / (freq1 * count1);
        //show(matrix_normal, n);

        //cout << "问题规模:" << n << endl;
        cout << "avx" << time1 << "ms";
        cout << endl;
        if (n > 320) size = 128;
    }

}

